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We investigate the transition from integrable to chaotic 
dynamics in the quantum mechanical wave functions from 
the point of view of the nodal structure by employing a two 
dimensional quartic oscillator. We find that the number of 
nodal domains is drastically reduced as the dynamics of the 
system changes from integrable to nonintegrable, and then 
gradually increases as the system becomes chaotic. The num- 
ber of nodal intersections with the classical boundary as a 
function of the level number shows a characteristic depen- 
dence on the dynamics of the system, too. We also calculate 
the area distribution of nodal domains and study the emer- 
gence of the power law behavior with the Fisher exponent in 
the chaotic limit. 

PACS number: 05.45. Mt, 64.60.Ak 



I. INTRODUCTION 

Quantum signatures of classically chaotic systems have 
been intensively studied and are now known to have some 
universal features in the energy level statistics. Simi- 
lar investigations on the signatures in the wave functions 
which may distinguish chaotic systems from integrable 
ones have also been performed [1,2]. A related signature 
of chaotic systems is given by the amplitude distribu- 
tion of wave functions which empirically reproduces the 
results of random matrix theory [3] . 

Recently, it was suggested in Ref. [4] that there is in 
fact such a universal character in the statistics of nodal 
domains of wave functions. The authors of Ref. [4] cal- 
culated the number of nodal domains Ni of the ith wave 
function, i.e., the regions where the wave function has a 
definite sign without crossing zeros, for two-dimensional 
billiards. They showed that the distribution of normal- 
ized number Ni ji of nodal domains for separable systems 
has a universal feature characterized by a square root 
singularity, while that for chaotic billiards shows a com- 
pletely different behavior, suggesting a scaling law NiK, i 
for large i. The latter scaling law has been derived in Ref. 
[5], where the authors adopted a percolation-like model 
to count the number of nodal domains. The authors of 
Ref. [5] analytically derived the scaling law for the aver- 
age number of nodal domains and showed that it agrees 
well with the numerical results for the superposition of 
random waves, i.e., a model which is supposed to sim- 
ulate the wave functions for chaotic billiards [1]. The 
percolation-like model allowed them to predict a power 



law behavior for the distribution of nodal domain areas, 
and also the fractal dimension of nodal domains. These 
predictions were shown to agree well with the numerical 
results for the superposition of random waves, although 
one may not conclude from these results alone that the 
area distribution provides a clear signature of quantum 
chaoticity. The number of nodal domains was studied 
also experimentally for the chaotic microwave billiard [6] . 

In the above studies of nodal domains, the two ex- 
tremes of dynamical systems, completely integrable (sep- 
arable) and chaotic (or its alternative), have been con- 
sidered. It is the purpose of the present paper to extend 
these studies to a more general nonintegrable system, 
where, by controlling a parameter in the Hamiltonian, 
one can interpolate the two extremes. We expect this 
would show the transition from integrable to chaotic sys- 
tems for the nodal domain distribution and thus may 
provide a clue to the role of nonintegrable perturbation 
which was implicit in the percolation-like model. We 
expect that a study of the power law behavior of nodal 
domain areas would also give a suggestion on the validity 
of the assumption adopted in the model. 

Below we first describe the model and the numerical 
procedure. We present numerical results for the distribu- 
tion of nodal domain numbers in Sec. Ill together with 
some analytical considerations. The distribution of the 
nodal intersections with the boundary of the classically 
allowed region is presented in Sec. IV. Results for the 
distribution of nodal domain areas are given in Sec. V. 
Proofs of formulas and some details of the calculations 
in the text are given in the Appendices. 



II. MODEL AND NUMERICAL PROCEDURE 

As a model which incorporates integrable as well as al- 
most chaotic systems, we adopt a two-dimensional quar- 
tic oscillator 
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H=^{pi + apl) + V{x,y), 

V{x,y) = }-{x^ + y^)-kx^y\ 



(1) 



where the parameter k controls the dynamics of the sys- 
tem. Detailed studies performed at a = 1 shows that 
the classical dynamics of the system at fc = 0.0 is inte- 
grable, becomes irregular as the value of k increases, and 
reaches an almost chaotic system at k — 0.6 [7]. The 
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energy level statistics of the quantum mechanical system 
show a similar transition, e.g., from Poisson to Wigner 
level spacing distribution as k increases [8] . The parame- 
ter a is introduced to break the symmetry with respect to 
the exchange of the x and y coordinates, which otherwise 
leads to an ambiguity in the definition of the eigcnfunc- 
tion at A; = 0.0. The value of a is set to the value 1.01 
throughout. We plot in Fig. 1 the Poincare surface of 
section for the system with a = 1.01 at several values of 
k. Qualitative behavior is almost the same as that with 
a = 1 
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FIG. 1. Poincare surface of section for the system of the 
Hamiltonian Eq. (1) with a = 1.01 at E = 5.0, x = 0.0, 
and (a) k = 0.0, (b) 0.2, (c) 0.4, (d) 0.6. The abscissa axis 
represents y coordinate and the ordinate py . 

In the numerical calculation, we first diagonalizc the 
Hamiltonian in a large harmonic oscillator (h.o.) basis 
and obtain wave functions for eigenstates. To deter- 
mine nodal domains we cut the two-dimensional sheet 
into small squares (meshes) so that they cover the al- 
lowed region of classical motion for each eigenvalue. We 
then study the sign of the wave function at each center 
of a mesh, and calculate the number of nodal domains 
by means of the Hoshen-Kopelman algorithm [9]. Here, 
two nearest neighbor meshes with the same sign arc con- 
sidered to belong to the same nodal domain. The size 
of the mesh is then changed to a smaller value until the 
convergence of the number of nodal domains is obtained. 
Adopted value of the mesh size is Xc\{Ei)/ max(1.4i, 200) 
for the ith eigenstate with energy E'j, where Xci{E) repre- 
sents the largest value of the x coordinate for a classical 
motion with energy E. One should note that the method 
becomes inaccurate for very small values of k, where the 
nodal crossing changes to a small avoided crossing due 
to the nonseparable perturbation. The smallest positive 
value of k in the present paper is 0.06 which was used in 
Sec. Ill C to study the qualitative behavior of the num- 
ber of nodal domains by comparing with a perturbative 



argument. 

Contrary to the case of the billiards, the size of the 
meshed sheet in the present case is not a priori deter- 
mined. In fact, since there is no hard wall, wave func- 
tions extend to infinity. They are, however, rapidly at- 
tenuated beyond classically allowed region for given en- 
ergy. Moreover, as shown in Appendix A, there appears 
no new nodal domain in the classically forbidden region. 
Therefore, we adopted the meshed sheet whose boundary 
coincides with that of the classically allowed region. All 
nodal domains, then, have an overlap with the meshed 
sheet. In case that the boundary of the meshed sheet 
cuts a nodal domain into several pieces the number of 
nodal domains may be overestimated due to the choice 
of the meshed sheet. To estimate the error caused by this 
boundary effect, we evaluated the difference between the 
number of nodal domains in the adopted meshed sheet 
and that in the square sheet that circumscribes the clas- 
sically allowed region. The estimated difference was 3.1% 
for k =0.1 and 2.2% for k =0.6. 

The Hamiltonian is still symmetric with respect to 
the a;-axis and the y-axis [7]. We calculated only those 
wave functions which are symmetric with respect to these 
two axes. The diagonalization space was truncated at 
fix + Uy < 200, where and Uy denote the numbers of 
oscillator quanta in the x- and y-directions. The adopted 
oscillator frequency in the diagonalization was optimized 
so as to minimize the value of TiH in this space for each 
k value [10]. 

When k = 0.0, we can obtain more accurate numeri- 
cal res\.ilts by adopting product of the wave functions for 
one-dimensional quartic oscillator calculated in a similar 
diagonalization procedure. Comparison of the number of 
nodal domains obtained in the two methods may provide 
a rough estimate for the accuracy of the adopted diag- 
onalization procedure, which in the present case is less 
than 2.4%. For the k = 0.0 results presented below we 
adopt those obtained by the product wave functions. 

III. DISTRIBUTION OF THE NUMBER OF 
NODAL DOMAINS 

A. General consideration on the number of nodal 
domains 

Before discussing the nodal structure for the present 

model, it is useful to study general properties of the num- 
ber of nodal domains which provide a useful classification 
scheme to be used in later sections. 

The number of nodal domains for a given nodal line 
structure of a wave function in a two-dimensional area 
with a boundary B is given by the formula 

N{nb,nc,m) = ^nb + ric + m + l, (2) 

under the assumption that more than two nodal lines 
never cross at the same crossing point. In Eq. (2), rib 
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represents the number of intersections of nodal lines with 
the boundary B, ric the number of crossing points of 
nodal lines, and m the number of islands. Here, the 
term 'island' means a cluster of mutually connected nodal 
lines which is linked neither to B nor to other clusters. 
A similar formula to (2) has been used and proved in 
Ref. [11] in a slightly different setting of the problem 
to study the multiplicity of eigenvalues for a membrane. 
In Appendix B we give a proof of the relation (2) for 
completeness. 

We discuss two cases in which this formula is especially 

useful. The first is the m = case, i.e., the case where all 
nodal lines are linked to the boundary, which is typical 
for separable (and therefore integrable) system. In this 
case 

N = ^rib + ric + 1. (separable, generic) (3) 

This allows us to discuss the dependence of N on the level 
number in relation to the structure of the wave function 
for A; = 0.0 as shown below. 

The second case is ric = which corresponds to the 
generic wave function of nonintegrable systems, where al- 
most all crossings of nodal lines change to avoided cross- 
ings except for an accidental case. We may first rewrite 
the number of nodal domains in Eq. (2) as the sum of 
two terms: 

AT = A/''" + A^^ (4) 

where A^'" denotes the number of 'inner nodal domains', 
i.e., those which do not touch the boundary B, while A^** 
denotes the remaining part, i.e., the number of ' bound- 
ary nodal domains' which touch the boundary B. If there 
is no crossing, ric = 0, it is easy to see that A^"^ is equal 
to the number m of islands. Prom Eq. (2), then, we find 
for ric = 

N™ = m, = ^rib + 1. (nonintegrable, generic) 

(5) 

The generic relation Eq. (5) for the number of boundary 
nodal domains may provide an estimate on the number 

of false crossing of nodal lines due to the finite mesh size 
in the numerical calculation. Comparison of the value of 
with the one obtained from the value of gives a 
difl[erence of 5.2% for k =0.1 and 2.5% for k =0.6. 
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FIG. 2. Number of nodal domains A'^; vs. the level number 
i for (a) k = 0.0, (b) 0.1, (c) 0.2, (d) 0.3, (e) 0.4, (f) 0.5, and 
(g) 0.6. Solid lines show the function f{i), while the dashed 
lines the function /(i). See text for details. 

B. Numerical results for the distribution of the 
number of nodal domains 

Figure 2 shows the distribution of the number of nodal 
domains, Ni, where the subscript i stands for the level 
number ordered according to the eigenvalue. Note that 
the quoted level number i is not the one of the total sys- 
tem, since we consider only levels which are symmetric 
with respect to the x axis and the y axis. In the present 
model, there are four symmetry classes, and the eigen- 
values in the four classes are almost equally distributed. 
In the space + Uy < 200 the total number of levels is 
20301 and that in the adopted symmetry class is 5151. 
Thus, we may assign for each i the corrected level number 
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i' approximately given by i' = ci with c = 20301/5151. 
The soUd Unes in Fig. 2 represent the prediction accord- 
ing to the percolation-like model which is given by the 
following function: 
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3^3-5 



(6) 



where the value of the coefficient /3 is taken as 2/7r as cal- 
culated for the billiard model [5]. In the percolation- like 
model, one first considers a nodal structure of the rectan- 
gular lattice pattern, and then assumes that an avoided 
crossing of nodal lines occur at every lattice point ran- 
domly, i.e., cither of the positive or negative domains is 
connected at the point with probability 1/2 independent 
of the other lattice points. Figure 3 shows a nodal do- 
main distribution represented as a histogram vs. Ni/i' . 
We give in Table I the average and the standard deviation 
of Ni/i'. 
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We first discuss the k = 0.0 case where the system 
is separable. In this case the eigenfunction is given by 
'>Pmn{x,y) = 4>m{x)(t>n{y), where (t>ni{x) represents the 
mth wave function of a one-dimensional quartic oscilla- 
tor in the x-direction. The behavior of the histogram in 
Fig. 3 (a) is similar to the one for the integrable sys- 
tem given in Ref. [4], showing an increase and a sharp 
cut-off at some value of Ni/i' . By inspecting Fig. 2 
(a) we find that this behavior comes from a number of 
regular sequences of eigenstates. The sequence with the 
largest Ni values is proportional to i, and the correspond- 
ing wave functions have the form ijjnn{x, y). This may be 
contrasted to the sequence with the smallest Ni values 
which is proportional to ^/i as shown by the dashed line 
in Fig. 2 (a). The wave functions in this sequence are 
of the form 4'ni{x,y) or tpin{x,y). The slopes of other 
sequences are intermediate between them. Typical nodal 
structures corresponding to tpnn{x,y) and tpni{x,y) are 
shown in Fig. 4. 
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FIG. 4. Nodal structures of wave functions at = 0.0: (a) 

An example of ipnn{x,y) {i = 98, = 169, = 56) and 
(b) an example of tp„i{x,y) {i = 103, A^]" = 0, Nf = 25). 



TABLE I. Average and standard deviation of Ni /i' distri- 
bution (total) and those of Nf^ /i' distribution (inner, see the 
text) for each k. Levels 1 < i < 1000 are considered. 

total inner 
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FIG. 3. Histogram oiNiji' for (a) k = 0.0, (b) 0.1, (c) 0.2, 
(d) 0.3, (e) 0.4, (f) 0.5, and (g) 0.6, where i' stands for the 
corrected level number. 
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FIG. 5. Number of inner nodal domains A']" vs. the level 
number i for (a) k = 0.0, (b) 0.1, (c) 0.2, (d) 0.3, (e) 0.4, (f) 
0.5, and (g) 0.6. Solid lines show the function f(i). 

These behaviors may be understood from Eq. (3). For 
eigenf unctions of the form ipni{x,y) there is no crossing, 



= 0, thus Ni 



/2. (Bclowwe explicitly write the 



i dependence of rife.) Since n^*'' is proportional to \/i [4], 
the number of nodal domains Ni in the lowest sequence is 
also proportional to The dashed line in Fig. 2 shows 
the function 



f(i) = dVi', d 



V3r(|; 



1.25, 



(7) 



which represents the number of nodal domains for wave 
functions il>ni{x,y) evaluated in a semiclassical way as 
given in Appendix CI. On the other hand, for eigen- 
functions of the type ipnn{x^y), the number of crossing 
is given by ric^ = (n^'')^/4, and Ni is dominated by nc'^ 
if n^^ is large enough. Accordingly, Ni in the highest 



sequence is proportional to i. Note that when the level 
number i becomes larger, in almost all levels the contri- 
bution of Uc becomes dominant, and the number of nodal 
domains A''^ becomes proportional to i on average. 

When the value of k becomes nonzero, e.g., fc = 0.1 
in Fig. 2 (b), the number of nodal domains is drastically 
reduced. This is due to the transition from the crossing 
to the avoided crossing of nodal lines. The shape of the 
histogram in Fig. 3 (b) changes, too: The peak at the 
highest value of Ni/i' at = 0.0 disappears and the 
histogram is now largely shifted towards low values of 
Ni/i'. 
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FIG. 6. Number of boundary nodal domains A',' 



the 



level number i for (a) k = 0.0, (b) 0.1, (c) 0.2, (d) 0.3,_(e) 0.4, 
(f) 0.5, and (g) 0.6. Dashed lines show the function f{i). 

To see more details, we plot A^i" and N^ in Figs. 5 and 
6, respectively. From these figures, we sec that at fc = 0.1 
the number of inner nodal domains drastically decreases 
below the f{i) line, while the distribution of the num- 
ber A^ of boundary nodal domains tends to accumulate 
around the f{i) line. These together lead to the reduction 
of the total number of nodal domains compared with the 
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separable case at fc = 0.0. The /c-dependence of iV'" and 
that of iVf are, however, quite different. As k becomes 
larger, the muiiber of boundary nodal domains decreases 
further and eventually becomes rather a small fraction of 
the total number of domains except at low energy (small 
i) region. In contrast, the number of inner nodal domains 
increases again, and in the chaotic limit around k = 0.6 
almost aligns to the f{i) line as in the case of the billiard 
system. Typical nodal structure for a wave function at 
k = 0.6 is shown in Fig. 7. 
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FIG. 7. Nodal structure of the 102nd wave function at 
k = 0.6, where A^]" = 30 and = 17. 

Let lis further study the behavior of the number of in- 
ner nodal domains N"^. We show the histogram of Nf^/i' 
in Fig. 8, and the average and the standard deviation of 
N™ /i' in Table I. The peak position of the histogram 
suddenly drops almost to zero as k becomes nonzero, 
which suggests that the nonintegrability first acts in such 
a way to eliminate inner nodal domains. This behavior is 
studied in the next subsection using perturbative argu- 
ment. The peak of N™ /i' then gradually increases with 
k, and finally the shape becomes approximate Gaussian 
at k = 0.6 centered at finite value of N™/i' , but still 
is much smaller than the case of fc = 0.0. Accordingly, 
the average of N™/i' increases. The behavior of the his- 
togram of total nodal domains at finite k in Fig. 3 follows 
that of the inner nodal domains. 

The number of boundary nodal domains iV^^ in Fig. 
6 also shows interesting features. At fc = 0.0 there are 
two regular sequences. The sequence with larger values 
corresponds to the wave functions Vmn = 4'm{x)4>n{y) 
(m 7^ l, n 7^ 1), while the sequence with smaller values 
to the wave functions ipni or ipin- In the latter wave 
functions nodal domains have a one-dimensional struc- 
ture, while in the former the boundary nodal domains 
are systematically aligned along the boundary. Both se- 
quences are proportional to the square root of the level 
number In Fig. 6, we show the fimction /(?) in Eq. 
(7) by the dashed hne. For small fc values, say k = 0.1, 
the number of boundary nodal domains is distributed 
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FIG. 8. Histogram of Nr/i'for (a) k = 0.0, (b) 0.1, (c) 0.2, 
(d) 0.3, (c) 0.4, (f) 0.5, and (g) 0.6, where i' stands for the 
corrected level number. 



C. Reduction of the number of nodal domains in the 
perturbative regime 

When the dynamics changes from integrable to nonin- 
tegrable, crossing of nodal lines changes to avoided cross- 
ing, which leads to the reduction of the number of nodal 
domains. However, this alone is not sufhcient to explain 
the difference in the reduction at small k value and that 
at large fc value as seen in Figs. 5 and 6. 

Figure 9 shows distributions of iV,, A^™, Nf, and the 
histogram of N™ /i' at k = 0.06. We see that the num- 
ber of inner nodal domains is smaller than the case of 
fc = 0.1, while the number of boundary nodal domains 
remains approximately the same. This further confirms 
that the number of inner nodal domains becomes smaller 
as the value of fc decreases as long as the value is nonzero. 
As noted above the avoided crossing becomes weaker as 
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the value of k decreases. Accordingly, the smaller mesh- 
sizc is required to avoid the fictitious crossing. Indeed, 
we used the mesh size Xc\{Ei)/ max(3.02, 200) for the cal- 
culation at k = 0.06. Properties of the avoided crossing 
(avoidance range) have been studied in Ref. [12]. 




FIG. 9. Number of different kinds of nodal domains at 
k = 0.06 is shown against the level number i : (a) Number of 
total nodal domains A'^;, (b) number of inner nodal domains 
Nl'^ , and (c) number of boundary nodal domains N^. His- 
togram of Nl'^/i' is shown in (d). Solid lines in (a) and (b) 
show the function f{i), while dashed line in (a) and (c) the 
function f{i). 

When the value of k is small, wave functions may well 
be approximated by the 1st order perturbation theory. 
The perturbed wave function can be written as, 

i^'mni^^y) = '^mn{x,y) + k ^ Cm'n''^m'n' {x,y) , 

{m' ■,n')^{m,n) 

{tpm'n'\x^y^\tpmn) 



(8) 



where Emn represents an unperturbed energy. Because 
of the second term, the crossing of nodal lines in the un- 
perturbed wave function ipmn will in generic case turn 
to an avoided crossing. Whether the positive domains 
or the negative ones are connected and merged into one 
nodal domain at this avoided crossing depends on the 
sign of the second term of (8) at the nodal crossing point 
of ipmn- Note that if this merging occurs randomly inde- 
pendent of the crossing points as in the percolation-like 
model, one would not obtain such a huge reduction of 
the number of inner nodal domains as seen in Fig. 9. 
This suggests that there is a correlation in the sign of 
the second term among different nodal crossing points of 
V'mn- The sign of the second term is mainly governed by 
the sign of the sum of a few principal components which 
have the largest values of [C^'n'!- For instance, we ver- 
ified that the probability of the sign of the sum of two 



principal components to be equal to that of the total sec- 
ond term at randomly chosen points (a;, y) is 0.86 for the 
present model. Moreover, the difference between the x 
quantum number of a principal component m' and that 
of the unperturbed wave function m, is generally very 
small compared with m: \5m\ = \m' — m| ^ m, as long 
as m is large. The same argument holds for the quantum 
number n in the y direction. 

We now evaluate the correlation between signs of a 
principal component with the x quantum number m' at 
two neighboring crossing points along the a;-direction. 
The distance d between two neighboring crossing points 
is the typical half wavelength for the quartic oscilla- 
tor wave function (pm- Suppose first that the distance d 
is shorter than the half wavelength ^A„i' of the principal 
component. Then, two neighboring crossing points must 
belong to the same half wavelength region in order that 
signs of the principal component at these two neighboring 
crossing points are same. Thus, the probability Pg that 
the signs of the principal components at two neighboring 
crossing points to be same can be estimated as 



P. = 



Am' A*] 



\dm\ 



An 



(9) 



where we used an approximate relation X„i oc m~^. 
When the distance d is larger than the half wavelength 
^Xm', the probability Pg can also approximately be given 
by Eq. (9). A superposition of a few principal compo- 
nents also has a wavelength close to the unperturbed one 
Xm, leading to the same result. Above argument can also 
be applied to the y-direction. 




FIG. 10. Nodal structure of the 92nd wave function at 
k = 0.06, where = and Nf = 27. 

Eq. (9) shows that the probability Pg is much less than 
unity. This means that nodal domains in the unper- 
turbed wave function tend to be connected along the di- 
agonal when a small perturbation is added, which results 
in the reduction of the number of inner nodal domains. 
Figure 10 shows a typical nodal pattern a.t k = 0.06. We 
see that unperturbed nodal domains, which are slightly 
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deformed from rectangular, are indeed connected along 
a diagonal direction. (Compare with Fig. 4 (a), where 
dominates.) This result supports the above consid- 
eration. 

For large k values, the perturbation theory becomes no 
longer a good approximation, and the above discussion 
cannot be applied. 



D. Model for the distribution of the number of 
nodal domains in the chaotic regime 

As a reference to the number of nodal domains, we used 
the function f{i) Eq. (6) derived for the chaotic billiard 
model [5]. As previously noted, in the percolation- like 
model, one first considers nodal structure of the rectan- 
gular lattice pattern. In order to determine the size of 
lattice, the relation = ky = E was used in Ref. [5], 
where kx denotes the wave number of a;-axis. This rela- 
tion is, however, specific to billiard models. Thus, if a 
system has a potential as in the present case, the distri- 
bution of the number of nodal domains may differ from 
billiard models even in the chaotic regime. Therefore, we 
would like to present a semiclassical model which pro- 
vides a more general expression for the nodal domain 
distribution incorporating some features of a Hamilto- 
nian with a local potential. 

Let us first summarize the results obtained in Ref. [4] 
for separable systems, where the Hamiltonian H is speci- 
fied by two action variables /i , /2 • According to them, the 
distribution function of Ni / i in a region Ei < Ei < E2 
(which we call A) in the I\ — I2 plane can be written as, 



1 



dhdhS ^ 



N{E) 



where N{E) is the semiclassical number of levels 



N{E) 



I 



dhdhe{E-H{h,h)), 



(10) 



(11) 



up to energy E, Na = N{E2) — N{Ei) the number 
of levels in the region A, and z^(/i,/2) = I1I2 is the 
number of nodal domains of the eigenstate specified by 
7i, /2- Assuming that the Hamiltonian is a homogeneous 
function of I\ and I2, we transform the variables from 
Ii and I2 to the energy E and r, where r represents 
the length of the path along the arc F with E = 1 
measured from the edge on the I2 axis. The action 
variables are written in terms of the new variables as 
(/i,/2) = {g{E)I^^\r),g{E)I^^\r)), where lf\r) and 

I^\r) represent the action on the arc F. When the 
Hamiltonian is a function of the ^th order in Ii and /2, 
g{E) = E'^^^. The Jacobian of the transformation is 
given by the product J{E)L{r), 



J{E) = 2g{E) 



dgjE) 
dE ' 



(12) 



L{r) = - 



dr 



dr 



(13) 



With these variables, the number of levels N{E) is given 
by 

N{E) = {g{E)f N{1), N{1) = j L{r)dr. (14) 

By performing the iJ-integration, one obtains the result 
of Ref. [4]: 



N{1) 



■ (15) 



In order to treat the nonintegrable cases, we modify the 
expression (10) so as to include the following features: 

i) Hamiltonian is not simply a function of 7i , I2 alone 
but depends also on the angle variables. 

ii) Number of nodal domains for a given energy is re- 
duced on the average due to nonintegrability. 

The first feature implies that each state is not specified by 
a point Ii , I2 in the phase space, but is distributed over 
an area with fixed E when projected on the Ii — 12 plane. 
We may include this effect by replacing the variable r 
involved in u with the smoothed one ra,{r); 



.(r) 



Jo 



Ur,r'ydr', 



(16) 



where, fu{r, r') represents a smoothing function with the 
width w, and rmax the total length of the arc F. The 
width (jo would correspond to the degree of the amplitude 
mixing of the wave function when expanded in the basis 
of the integrable system. For the integrable system, w = 
0, and when the system is chaotic, the wave function 
will be distributed over the available phase space, i.e., 
w > ''max- The second feature ii) may be included by an 
introduction of a reduction factor G, i.e., 



(17) 



with G = 1 for separable systems. Thus in our model, 
the distribution is given by 



/f)(r„(r))4°)(r„(r)) 



N{1) 



X G{l't\r^{r))g{E),l!^\r^{r))g{E))) . (18) 

To perform the integration, we have to specify G. In 
the chaotic limit, the asymptotic value of G with large 
Ii and I2 values may be obtained from the percolation 
model of Ref. [5] as 



8 



oo. 



(19) 



On the other hand, uj becomes large in the chaotic Hmit. 
The function wih then become independent of r and 
is represented by the average, i.e., f = rinax/2. By sub- 
stituting this into Eq. (18), we obtain 



p{0 = H^-^^ ^ 



N{1) 



(20) 



where ^ corresponds to the average value of Ni/i. 

We may adopt a concrete example to obtain the value 
of ^. Let us consider as an unperturbed model Hamilto- 
nian the following form: 



Ho oc /f + 4. 
The average is calculated as. 



■GP, 



where we used (r = rmax/2) 



4°)(r-)=7(°)(r-) = 



and 



iV(l)= [\l-x^)Ux= ^^^^ 
Jo 



(21) 



(22) 



(23) 



(24) 



Eq. (22) can be applied to the Hamiltonian with or with- 
out a potential. For examples, the average for the billiard 
model {£ = 2) is 



3^/3-5 



0.062, 



(25) 



and that for the quartic oscillator model (^ = 4/3) is 



2_7r(3\/3_-5) 

6r(!)' 



0.055. 



(26) 



The difference between the average for the billiard model 
and that for the quartic oscillator model is quite small, 
and cannot be recognized in the log-log plot like Fig. 2. 
Figure 11 compares in detail the result of the numeri- 
cal calculation and the above two asymptotic values. In 
the numerical calculation we employ the number of inner 
nodal domains N™ instead of Ni in order to remove the 
boundary effects. The average by numerical calculation 
still fluctuates between two asymptotic values and does 
not seem to converge. We must investigate higher levels 
to see if the actual average converges to the value Eq. 
(26). 

In the intermediate region of integrable and chaotic 
cases, the results for the distribution of nodal domain 



numbers depends on the concrete form of the weight- 
ing function and the reduction factor G. Here, we 
consider the specific case where the amplitude mixing is 
not complete although large, while the correlation among 
avoided crossings is lost. To be more definite, we take the 



range of uj as 



< UJ < r„ 



and assume that the 



reduction factor G can be obtained by the percolation 
model. In Appendix D, we show that under a reasonable 
choice of f^, we obtain the result: 



d 



> 0. 



(27) 



In order to find which values of k correspond to such a 
case, we have to know the relation between k, w, and G, 
which remains for a future study. However, the numer- 
ical results shown in Sec. V imply that for k > 0.4, the 
assumption of the percolation-like model may be applied. 
Eq. (27) suggests that the average value of Ni/i increases 
as one approaches the chaotic system, which is in accord 
with the results of Table I. 
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FIG. 11. The average {Nl^^/i'} for k = 0.6. Averaging is 
performed for bins of i each of which contains 50 levels, and 
the solid squares are placed at the center in each bin. The 
solid line shows the asymptotic value for the quartic oscillator 
model and the dashed line for the billiard model. 



IV. DISTRIBUTION OF THE NUMBER OF 
NODAL INTERSECTIONS 

We now consider the behavior of the number of bound- 
ary nodal domains, in particular, its dependence on the 
level number ?', which can be obtained from that of the 

(i) 

number of nodal intersections with the boundary, nl , 
according to Eq. (5). 

For the case of billiard systems, the number of inter- 
sections is proportional to Vi, even if the dynamics of 
the system is chaotic [4]. On the other hand, in Ref. [13], 
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the behavior of the nodal intersections with the bound- 
ary of the classicahy aUowed region was studied for sys- 
tems with soft potentials whose dynamics is chaotic. Ac- 
cording to Ref. [13], the number of intersections per unit 
length along the boundary is given by 

p = 7|gl■ady|l/^ (28) 

where 7 = 0.171 and V represents the potential. 




200 400 . 600 800 1000 

I 



FIG. 12. (a) Average of n^'^ for k = 0.0, 0.1, 0.2, 0.3, 
0.4, and 0.5 from top to bottom. Averaging is performed as 
in the previous figure. Solid and dashed lines correspond to 
r(i) ~ 5.6i^^^ and 5.6i^^^, respectively, (b) Average of n['^ 
for k = 0.6. The dashed line corresponds to 5.6i^^^. 

Now we calculate n^^^ for the present model, Eq. (1), 
by integrating p along the boundary V{x,y) = Ei. As 

(i) 

the potential is homogeneous, can be written as. 



= 7/ ds\gia.dV\^/^ = jbEj 

Jv=Ei 



'V=Ei 

with the coefBcient b given by 



1/2 

i ' 



b= (f rfs|grady|^/^ 

Jv=i 



(29) 



(30) 



where the integral is performed along the curve V{x,y) = 
1. By using the semiclassical relation between the energy 

Ei and the level number i given in Eq. (14), we obtain the 
level number dependence of the number of intersections: 



W _ ,-/l/3 _ o r) -/l/3 _ 4 o -1/3 

~ 7V(1)V3* ~ ' 



(31) 



where i' represents the corrected level number and the 
numerical value b ~ 15.26 evaluated at A: = 0.6 was used. 

We note that nj,'-' has a Vz dependence at fc = 0.0. 
This indicates that the level number dependence of the 
number of intersections changes from i^/^ to i^/^ as the 
dynamics changes from integrable to chaotic. Numerical 
results shown in Fig. 12 confirm this. Here, the average 
value of n^*' (over 50 levels) is shown as a function of 

the level number i. At k = 0.0 n^'^ follows well the 

semiclassical result for the number of intersections shown 
by the higher curve in Fig. 12(a), i.e.. 



Tii) = 



8^5/4 



5V3r(2^2T-^5 



r(|) 



Vi' ~ 2.8Vi' = 5.6\/f, (32) 



a derivation of which is presented in Appendix C 2. As 

(i) 

k increases, gradually decreases. Finally, at A: = 0.6, 



n['^ follows the curve 5.6i^/^ shown as the dashed line in 
Fig. 12 (b), where the fitted value is used as a coefficient. 
The fitted value 5.6 is a little larger than the predicted 
value 4.8. This discrepancy may be due to the curvature 
of the boundary V{x,y) = 1, which is not taken into 
account in Eq.(28). 



V. DISTRIBUTION OF THE AREA OF THE 
NODAL DOMAINS 

Let us now consider the distribution of the area of 
nodal domains. The area distribution may provide a 
measure as to whether the percolation-like model can be 
applied. The appropriate length unit to define the area 
should be the wave length [5] which depends on the en- 
ergy of the level. Moreover, in contrast to the billiard 
problem, the wavelength in the present model depend 
locally on the coordinate due to the presence of the po- 
tential V{x,y). We define the scaled area s of the nodal 
domain for the level i with energy Ei by 



■/ 

Jn.d. 



s = n I {Ei- V{x, y))dxdy, 



(33) 



where the integral is performed over each region of the 
nodal domain. We take the scale factor = 1. One may 
note that the parameter a is absent in the definition (33) 
so that it does not directly influence the properties of 
the area distribution. The normalized number of nodal 
domains with area s is defined as [5], 



n{s) ^ Yl 



Q^is) 



(34) 



Here, Qi{s) represents the number of nodal domains with 
area s in the ith wave function. 

Figure 13 shows the distribution of n{s). The line 
in the figure represents the power distribution with the 
Fisher exponent, n^(s) oc s~'^, where r = 187/91, which 
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wc call hereafter the Fisher line. This represents the 
characteristic distribution at the critical point in the two- 
dimensional percolation model. The figure shows that 
the distribution at small values of k deviates consider- 
ably from the Fisher line, especially in the small area 
region. The deviation gradually diminishes as the value 
of k increases, and at fc = 0.6, the distribution almost 
coincides with the Fisher line except for a small discrep- 
ancy in the small s region. We now study if the latter 
discrepancy may imply that the assumption behind the 
percolation-like model does not hold completely even at 
k = 0.6. 
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Figure 14 shows the area distribution n{s) for inner 
nodal domains. Large deviation from the Fisher line at 
small k values reflects the fact that the avoided crossings 
are correlated and the nodal structure does not belong to 
the universality class of the percolation model. On the 
other hand, when k > 0.4, the agreement is very good. 
Thus we may conclude that the small deviation seen in 
Fig. 13 at A; > 0.4 is due to the boundary effect. 
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FIG. 14. Same as Fig. 13 but for the inner nodal domains. 
The Fisher line corresponding to fc = 0.0 is omitted. 



FIG. 13. Distribution n{s) for k = 0.0, 0.1, 0.2, 0.3, 0.4, 

0.5, and 0.6 from top to bottom. The values of n(s) arc 
muhiplied by a factor 10^^ for k = 0.0, lO"' for k = 0.1, 10** 
for k = 0.2, 10® for k = 0.3, lO"* for k = 0.4, and 10^ for 
k = 0.5. Levels 200 < i < 1000 are considered. Solid lines 
show the Fisher line. 

We note that the Fisher line was obtained for an in- 
finite lattice of mesh points, while in our finite system 
the effect of boundary may influence considerably the 
area distributions. In order to compare with the infinite 
system and to test the applicability of the percolation 
model, we have to exclude the effect of the boundary. To 
this end, it is appropriate to employ only the inner nodal 
domains to obtain the area distribution. 



Remember that the distribution of the number of nodal 
domains in Fig. 3 and the number of inner nodal do- 
mains in Fig. 8 still show a gradual change after k = 0.4. 
The area distribution is contrasted to this behavior. The 
agreement with the Fisher line; alrc;ady at A; > 0.4 may 
imply that the independence and the randomness of the 
avoided crossings is practically already realized before the 
system becomes completely chaotic at fc = 0.6 as given by 
the classical phase space structure or by the level statis- 
tics. 



VI. SUMMARY 

We investigated the transition from integrable to 
chaotic dynamics from the point of view of the nodal 
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structure in the wave functions by employing a two di- 
mensional quartic oscillator. The distribution of the 
number of nodal domains, the number of intersections of 
nodal lines with the boundary of the classically allowed 
region, and the area distribution of nodal domains were 
studied as a function of the parameter k which controls 
the nonintegrability of the system. 

The number of nodal domains is drastically reduced 
as the dynamics of the system changes from intcgrable 
(fc = 0.0) to nonintegrable (fc nonzero), and then gradu- 
ally increases as the system becomes chaotic {k ~ 0.6). 
Separation into 'inner' and 'boundary' nodal domains 
shows that the above dependence on k mainly comes from 
the behavior of the former. Perturbative argument sug- 
gests that a finite (but small) k gives rise to correlated 
avoided crossings of nodal lines, leading to a drastic re- 
duction of inner nodal domains. Their number turns to 
increase again as k becomes larger, which may be related 
to the loss of the correlations in avoided crossings which 
underlies the percolation model. The 'boundary' nodal 
domains show, in contrast, a milder dependence on fc, 
and its significance in the total number of nodal domains 
becomes small at large k values. A scmiclassical model 
which incorporates the degree of amplitude mixing as 
well as properties of avoided crossings has been proposed 
to study the number of nodal domains. 

We studied the distribution of the number of intersec- 
tions with the boundary of the classically allowed region. 
We found that the average number shows a different de- 
pendence on the level number as the dynamics changes 
from integrablc to chaotic. It is interesting to find such 
a characteristic connection to the dynamics in the struc- 
ture of the wave function at the boundary, in view of 
the rather mild k dependence found in the number of 
boundary nodal domains. 

We studied also the distribution of the nodal domain 
areas which shows a scaling behavior in the percolation- 
like model [5]. In the present model, too, the area distri- 
bution shows a scaling with the same exponent for large 
k values. A small deviation has been shown to come from 
the boundary effect. The scaling behavior seems to be 
complete before the dynamics reaches the chaotic limit, 
however. 

This raises the question as to how the various signa- 
tures on the chaotic properties which appear in the level 
statistics or in the wave functions are interrelated. Fur- 
ther studies on the typical values of k which characterize 
the onset or the completion of various structures in the 
wave functions would be interesting. They include, for 
instance, the signals studied in the nodal domain distri- 
bution, amplitude mixing (Porter-Thomas distribution), 
loss of correlations in the avoided crossings, scaling in the 
area distribution, etc. These are left for the future inves- 
tigation. In this context, it is valuable to mention that 
even for the superposition of random waves, there is a 
weak long range correlation among the avoided crossings 
as discussed in Ref. [14]. 
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APPENDIX A: NO NODAL DOMAIN IN THE 
CLASSICALLY FORBIDDEN REGION 

We prove that there is no nodal domain which is em- 
bedded entirely in a classically forbidden region. 

Assume that there is such a nodal domain whose region 
is D for an eigenfunction ^E" (taken real, for simplicity). 
The eigenfunction 4" satisfies the following Schrodinger 
equation; 

l^^ = 2{y - E)<^ , (Al) 

where E represents the eigenvalue and V the potential. 
Multiply ^ on both sides in Eq. (Al) and integrate them 
over the region D. The right hand side becomes positive, 

2 [ {V- E)^'^dv > 0, (A2) 

JD 

because in the classically forbidden region the inequality 
V — E > Q always holds. On the other hand, since the 
value of ^ along the boundary of the region D is zero, 
the left hand side after partial integration becomes 

/ *A*dt; = - / (V*)^dw < 0, (A3) 
Jd Jd 

which is in contradiction to Eq. (A2). Therefore there is 
no nodal domain in the classically forbidden region. 

APPENDIX B: NUMBER OF NODAL DOMAINS 
FOR A GIVEN DIAGRAM OF NODAL LINES 

We derive the relation (2) for the number of nodal do- 
mains in a two-dimensional area enclosed by a boundary 
B. See also Ref. [11] for a slightly different form of the 
relation for the number of nodal domains expressed in 
terms of the number of nodal lines, etc. The nodal lines 
and the boundary represent a kind of a diagram in a two- 
dimensional plane which is constructed with a number of 
vertices (nodal crossings including nodal contacts at the 
boundary) that arc connected with edges, i.e., the line 
segments of nodal lines or those of the boundary. We 
first consider a general diagram made of nodal lines and 
then restrict it to the special case treated in the text. 

Let us define the degree A:(> 3) of a vertex as the num- 
ber of lines which are connected at the vertex. We define 



12 



also 'island' (in B) as the cluster of nodal lines which 
are linked neither to B nor to the other clusters of nodal 
lines. Simplest island is a 'bubble', i.e., a closed line with- 
out a vertex. For instance, a bubble within a bubble will 
make two islands. 

We now consider the number N of nodal domains for a 
diagram of nodal lines in B with m islands, where the to- 
tal number of vertices of degree k{> 3) is rifc. Since there 
are m + 1 clusters of nodal lines disconnected from each 
other, we assign them an index j = 0, 1, • • • , m, where 
j = denotes the cluster which involves the boundary 
B, while j = 1, • • • , m denote islands. For each cluster 
j we assign the number of nodal domains Nj, the total 
number of edges ej, and that of vertices Vj which is given 
by the sum of nk{j), the number of vertices with degree 
k in the j-th cluster. If we neglect the case of 'bubbles' 
for a moment, we can use Euler's formula for each cluster 
j to obtain 

Nj + 1 = ej - Vj +2, (,7 = 0, 1, 2, . . . , m) (Bl) 

where we added unity on the left hand side, since Nj 
counts domains only inside the boundary. By using the 
graphical relation 

^o^\Y.^Mj), (B2) 

fe>3 

the relation (Bl) is transformed to 

fe>3 ^ ^ 

Equation (B3) is seen to hold also for a 'bubble', where 
N = 1 and all Uk's are zero. By summing up over j we 
obtain the relation 

N = f2^j = Y.(l'^-'^)^k + m+l, (B4) 

j=0 k>3 ^ ^ 

where we used = 'rik{j)- 

In the special case treated in Sec. Ill, we assumed that 
there are no accidental crossing of more than two nodal 
lines at a point inside B, and also that no accidental 
contact of more than one nodal line on iJ. In this case, 
the vertex inside B has a degree 4, and the vertex on 
B has a degree 3. Thus the number of nodal domains 
N[b, c, m) for a diagram with c(= 77.4) crossings, b{= 713) 
points of contact on B, and m islands is given by 

N{b,c,ni) = + c + m + 1, (B5) 

which is the relation (2). 

APPENDIX C: DERIVATION OF SOME 
RELATIONS 

We give brief derivations for Eqs. (7) and (32). See 
Sec. HID for notations. 



1. A derivation of Eq. (7) 

For the wave function ipni {x, y) which has correspond- 
ing actions, Ii = I, I2 = 0, the semiclassical number of 
nodal domains Ni can be written as, 

Ni = I = g{E), (CI) 

while the semiclassical level number is given by Eq. (14). 
Then, the coefficient d in Eq. (7) is given by 

d = l/^/N{T), (C2) 

which leads to Eq. (7) for the quartic oscillator model. 

2. A derivation of Eq. (32) 

Analogous to Eq. (10), the distribution function of 
/Vi can be written as, 

where t(/i,/2) — 2(/i + 12) is the number of nodal inter- 
sections. By the same variable transformation as in Sec. 
Ill D and by performing the integration over the variable 
E, the average of 77 can be written as, 

('^> = J^(^J/'-^('-)(^^°^('-) + ^2\r)). (C4) 

Inserting Eq. (13) and transforming the variable r to 
Ii^\= x), the integral in Eq. (C4) can be expressed as, 

J dxx{l-x^)i+ J dx{l-x^)^ 

1 £(i)im + £(i)im 

ff(^) + ^r(^) ' ^^'^ 

for the Hamiltonian with the form Eq. (21). Thus, for 
the quartic oscillator model {£ = 4/3), Eq. (C4) gives the 
coefficient in Eq. (32). 



APPENDIX D: EXAMPLE FOR THE 
CALCULATION OF THE AVERAGE OF NODAL 
DOMAIN NUMBERS 

We here show the increase of average oi Ni/i after ui > 
^''max by using Eq. (18) with the assumption that the 
reduction factor G can be obtained from the percolation 
model. The reduction factor G{I\,l2) may still have a 
dependence on I\ and /2, however, when the energy is 
not so high. 

As for the smoothing function, we adopt the following 
form: 
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9{lj 



r'\) 



min(r + a;, rmax j — niax(r — uj, Oj 
The average of Ni/i can be written as, 



(Dl) 



of the reduction factor G will be close to the asymptotic 
one. Thus, the derivative of G is near zero; 



d_ 
da 



G(a,6) = ^G(a,6)^0. 



(D8) 



(D2) 



Accordingly, the second term in Eq. (D7) is neglected. 
We also assume the Hamiltonian for integrable case has 
the form, Eq. (21). The function G{E,r) ,then, can be 
approximately written as, 



where r^^ = ri^{r) and g = g{E). We differentiate Eq. 
(D2) with respect to w, 



C(E,r)^7r/f|l-|^ 



du \ i 



x^C(E,r) (D3) 

G{E,r) = (/j°^r.)/f (r.) +/f (r.)lf' (r.))G 

+4°)(r.)/f(r.),(^/r' + f/f'), (D4) 




(D9) 



Since in the considered range of r^j, the relation if'^ < 
I^'' always holds, the derivative of average with respect 
to ix> is positive. 



diij \ I 



(DIO) 



where r'^ = dr^{r)/duj, /f^'(r) = rf7}"^(r)/dr, and the smoothing width w', as long°as w > 
dG / da = dG(a,b)/ da etc. , . , , , . 

When UJ > ir„,ax, the function has a following form; i^^^^^^s^ ^""^'^^^ ^^ops when the condition ![ 



r(0) 



2 ' max • 



This result shows that the average of Ni/i increases as 

The 
(0) _ 



0, (rmax - w < r < w) 



12"^ always holds, which means r,^ = ^^max, namely the 
chaotic limit. 



2' 



(r > w). 



(D5) 



Accordingly, the integration range of r in Eq. (D3) is 
restricted to (0,r„iax — ^) and (t^^, rmax)- By using the 
relation G{a, b) = G{b, a) and approximate relations 



l['\r)=ir{r, 



(0), 



r) and l[^^'{r) 



-'2 \ max 



r), 



the latter integration range can be transformed to the 
former as, 



d_ 



1 



E2 



J{E)dE 



NaN{1) je, 

drL{T)G{E,r). 



(D6) 



Above approximate relations are justified by an approxi- 
mate symmetry between the variables l\ and I1 , because 
the symmetry breaking in the Hamiltonian (1) is small. 

The behavior of the average is, thus, determined by the 
function C{E,r). The function C{E,r) can be written 
as, 



C{E,r)=ir4'U 1 + 



I^dli'^' 

Since the range of is jr^ax < < ^r^ax 



(D7) 



for UJ > 
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